SARS-CoV-2 infection results in immune responses in the respiratory tract and peripheral blood that suggest mechanisms of disease severity

Respiratory tract infection with SARS-CoV-2 results in varying immunopathology underlying COVID-19. We examine cellular, humoral and cytokine responses covering 382 immune components in longitudinal blood and respiratory samples from hospitalized COVID-19 patients. SARS-CoV-2-specific IgM, IgG, IgA are detected in respiratory tract and blood, however, receptor-binding domain (RBD)-specific IgM and IgG seroconversion is enhanced in respiratory specimens. SARS-CoV-2 neutralization activity in respiratory samples correlates with RBD-specific IgM and IgG levels. Cytokines/chemokines vary between respiratory samples and plasma, indicating that inflammation should be assessed in respiratory specimens to understand immunopathology. IFN-α2 and IL-12p70 in endotracheal aspirate and neutralization in sputum negatively correlate with duration of hospital stay. Diverse immune subsets are detected in respiratory samples, dominated by neutrophils. Importantly, dexamethasone treatment does not affect humoral responses in blood of COVID-19 patients. Our study unveils differential immune responses between respiratory samples and blood, and shows how drug therapy affects immune responses during COVID-19.

S ymptoms of SARS-CoV-2 infection, known as coronavirus disease 2019 (COVID-19), vary from asymptomatic or mild disease to critical illness, including respiratory failure and death 1 . Global efforts focused on developing new drugs and vaccines. While vaccines showed safety and immunogenicity towards SARS-CoV-2 2-4 , the effects of drug treatments remain controversial. Dexamethasone, a synthetic glucocorticoid drug, can lower the 28-day mortality rate in COVID-19 patients receiving oxygen support, prolong ventilator-free days and improve oxygen partial pressure to fractional inspired oxygen (PaO 2 /FiO 2 ) ratios compared to placebo or standard care [5][6][7] . However, SARS-CoV-2 RNA can be detected for longer in patients receiving glucocorticoid treatment 8 . Treatment with remdesivir, a nucleoside analogue inhibiting RNA-dependent RNA polymerase (RdRp), in COVID-19 can shorten the time to recovery and provide better clinical outcomes [9][10][11] . However, the effects of dexamethasone and/or remdesivir on humoral and cellular immune responses are unclear.
To dissect the breadth of immune responses during SARS-CoV-2 infection in the respiratory tract compared to those detected in blood, we collected paired longitudinal blood and respiratory samples from hospitalised COVID-19 patients to investigate innate, adaptive and humoral immunity. Overall, our study unveils differences and defines correlations in innate and adaptive immune responses between respiratory and blood samples of COVID-19 patients and provides insights into potential biomarkers and immunotherapies for severe COVID-19.

COVID-19 patient cohort.
To define immune responses to SARS-CoV-2 in the respiratory tract, we obtained 41 respiratory samples (15 endotracheal aspirates (ETA; from 11 patients), 20 sputum samples (from 18 patients), 6 bronchoalveolar lavage samples (BAL; from 6 patients)). Respiratory samples were collected from 33 PCR-positive COVID-19 patients from whom we also collected 34 paired blood samples (Fig. 1a, Supplementary Table 1 and Supplementary Table 2). Three COVID-19 patients were admitted to the ward while 30 patients were in the ICU ( Fig. 1a; Supplementary Table 1). The median age of COVID-19 patients from whom we obtained respiratory samples was 55 years (range 25-76) and 33.3% were female (Supplementary Table 1). Where feasible, blood was collected on hospital admission, during hospital stay and on hospital discharge. No significant differences were found between time of respiratory specimen collection and matched blood samples collected at the closest time-point (p = 0.89; Fig. 1b).
To determine the effects of dexamethasone, an antiinflammatory corticosteroid, taken alone or in combination with the anti-viral drug remdesivir, on immune responses in blood, we recruited 57 COVID-19 patients (42 ward patients and 15 ICU patients) with a median age of 58 years (range 22-90) and 49.1% female from whom we obtained 86 blood samples (Fig. 1a, Supplementary Table 1, Supplementary Table 2).
SARS-CoV-2 genome sequence data, available from 40 out of 84 COVID-19 patients, showed that patients were infected with SARS-CoV-2 viruses belonging to the transmission network (TN)-A, representing a highly clonal and dominant network during the second wave of COVID-19 epidemic in Victoria 22 , except for 1 patient belonging to TN-B (#001) (Fig. 1d).
ICU admission associated with higher NIH severity score, oxygen therapy, drug treatment and weight. Disease severity within our cohort was stratified according to whether COVID-19 patients were hospitalized in the ward or ICU. COVID-19 patients were also graded according to the NIH severity score of 1-5 according to their symptoms (Supplementary Table 3). ICU patients had significantly higher NIH scores compared to ward patients (p < 0.0001, Fig. 1e) and more ICU patients received oxygen support (p < 0.0001) and drug treatments, either dexamethasone alone or dexamethasone with remdesivir (p < 0.0001) (Fig. 1e, Supplementary Table 2). Interestingly, ICU patients also had significantly increased body weight (p = 0.0008), height (p = 0.0476) and body mass index (BMI) (p = 0.0085). Age correlated with the length of hospital stay (p < 0.0001, Fig. 1e), but no differences in age, gender, ethnicity, immunosuppressant drugs or smoking were observed between ICU and ward patients (Fig. 1e, Supplementary Table 2).
Differential inflammatory cytokine profiles in respiratory and plasma samples. There are scarce data on the inflammatory milieu in COVID-19 respiratory specimens. To determine cytokine/chemokine levels and composition in respiratory samples compared to paired plasma, we measured cytokines/chemokines (IL-1β, IFN-α2,  IFN-γ, TNF, MCP-1, IL-6, IL-8, IL-10, IL-12p70, IL-17A, IL-18,  IL-23 and IL-33), sIL-6Rα and an extracellular matrix protein "disintegrin and metalloproteinase with thrombospondin motifs-4" (ADAMTS4) 23 . Amongst the COVID-19 patients, greatly elevated levels of inflammatory cytokines/chemokines were detected in respiratory samples across ETA, sputum and BAL specimens, with concentrations being 160× (MCP-1), 90× (IL-6) and 110× (IL-8) Fig. 1 Demographics of the COVID-19 cohort. a Number of patients recruited in the COVID-19 cohort and samples collected are shown. b Comparison between the time of respiratory sample (endotracheal tube aspirate (ETA), sputum, and bronchoalveolar lavage (BAL)) and paired blood sample collection are shown and were calculated using a two-sided Mann-Whitney test. n Respiratory = 40, n Blood = 33. c Collection time of COVID-19 blood samples. d Maximum likelihood phylogenetic tree of SARS-CoV-2 sequences from Victoria from 28 January 2020 to 28 October 2020 (including context sequences from the rest of Australia and New Zealand). Phylogenetic tree includes randomly subsampled sequences from transmission networks (TN) A and TN B in Victoria, with a total number of 10941 and 145 cases respectively. The outermost tip of each radial line represents a single sequence; the sum of each radial line between two tips represents the genetic distance between two sequences. Each radial stepwise progression represents approximately one single nucleotide polymorphism (SNP). Sequences from study patients (n = 40) are shown as open circles (patient with blood samples only) or solid-coloured circles (patients with blood and respiratory samples). Half-filled circles are used when samples are located close to each other. e Distribution of clinical data in ward and intensive care unit (ICU) COVID-19 patients. n Ward = 45, n ICU = 39. The patients received dexamethasone (D), dexamethasone with remdesivir (D + R) or neither (N). The bounds of the box plot indicate the 25th and 75th percentiles, the bar indicates medians, and the whiskers indicate minima and maxima. Statistical significance was determined with a two-sided Fisher's exact test for the National Institutes of Health (NIH) score, and a two-sided Mann-Whitney test for age, weight, height, and body weight index (BMI). Correlation was determined with a two-tailed Spearman's correlation. Source data are provided as a Source Data file.
Since the magnitude of cytokines/chemokines was much higher in respiratory samples than in plasma, a z-score normalization was performed for respiratory and matched plasma samples separately ( Fig. 2d). High normalized cytokine levels were detected in ETA samples from only 2 out of 11 COVID-19 patients (#26 and #49) (Fig. 2d). In sputum samples, however, this was seen in 10 out of 18 patients (i.e. higher than ETA and BAL samples; #21, #72, #73, #74, #75, #76, #80, #81, #82 and #84). Our data therefore suggest that sputum potentially represents the more desirable specimen type that reflects the high inflammatory milieu at the primary site of SARS-CoV-2 infection. Conversely, the majority of patients displayed elevated cytokine/chemokine levels within plasma samples. Overall, while the inflammatory cytokine/chemokine levels were excessively higher in respiratory fluid compared to plasma, they were variable across COVID-19 patients, indicating that the plasma inflammatory milieu does not always reflect the airway inflammation and that hospitalized/ICU COVID-19 patients should be monitored for inflammation in airways, such as in sputum, to understand disease severity and potential benefits of immunomodulatory treatments.
High RBD-specific IgM and IgG seroconversion in COVID-19 respiratory samples. SARS-CoV-2-specific antibodies in respiratory samples are relatively unexplored. We measured SARS-CoV-2 RBD-specific IgM, IgG and IgA antibodies in paired respiratory and blood samples using RBD-ELISA and surrogate virus neutralisation test (sVNT) (Fig. 3). Compared to non-COVID-19, COVID-19 patients displayed higher levels of RBD IgM (p = 0.0003) and IgG (p < 0.0001), but not IgA, in respiratory samples (Fig. 3a, b), which was possibly either due to technical issues or cross-reactivity of IgA antibodies (Fig. 3b top left panel). However, significantly lower titres of RBD IgM and IgG were found in COVID-19 respiratory samples compared to matched plasma samples (Fig. 3b top right panel). This was consistent for both pooled respiratory samples ( Fig. 3b top right panel) as well as separately analysed ETA, sputum and BAL samples (Fig. 3b bottom panel), with the exception of ETA IgG titres which had a lower median than matched plasma samples, though not significant.
Using sVNT, more sputum than ETA or BAL samples had detectable neutralizing activity, associated with high levels of RBD-specific IgM and IgG antibodies (Fig. 3c, d). Neutralizing activity was not detected in the majority (58.5%) of respiratory samples at the acute time-points. This included 11 ETA, 11 sputum and 2 BAL samples. Plasma samples with high neutralizing activity had high levels of all three Ig isotypes of RBD-specific antibodies, which positively correlated with neutralizing activity (Fig. 3c, d). This was also observed in respiratory samples, though only anti-RBD IgM and IgG significantly correlated with neutralizing activity (Fig. 3c, d). Seroconversion levels of RBD-specific IgM and IgG antibodies were detected in the majority of COVID-19 respiratory samples (34/41, 83%) and patients (26/33, 79%) (Fig. 3e), suggesting the prominence of RBD-specific IgM and IgG in respiratory samples during acute COVID-19. In terms of correlations between respiratory samples and plasma, overall IgM, IgG and sVNT levels correlated across the specimens (Fig. 3f).
To investigate the most prominent antibody features that differed between COVID-19 and non-COVID-19 patients, Partial Least-Squares Discriminant Analysis (PLSDA) was performed (Fig. 4c). As few as 3 antibody features were sufficient to separate COVID-19 and non-COVID-19 ETA, with higher SARS-CoV-2specific IgG and IgM in COVID-19 ETA, consistent with higher anti-RBD IgG and IgM in ELISAs (Figs. 3b, 4c). In contrast, higher SARS-CoV-2-S-specific IgG antibodies and antibodies Fig. 2 Discordant levels of cytokines and chemokines in COVID-19 respiratory samples compared to paired plasma samples. a Absolute concentrations of 13 cytokines and chemokines (IL-1β, IFN-α2, IFN-γ, TNF, MCP-1, IL-6, IL-8, IL-10, IL-12p70, IL-17A, IL-18, IL-23, and IL-33) in pooled respiratory and paired plasma samples. b Comparison of cytokine and chemokine levels between endotracheal tube aspirate (ETA), sputum or bronchoalveolar lavage (BAL) and paired plasma samples using a two-sided Mann-Whitney test. Bars indicate the median values. c Correlation of cytokine and chemokine levels between respiratory samples (ETA, sputum, and BAL) and paired plasma samples collected at the closest timepoint for each patient. Correlation was determined with a two-tailed Spearman's correlation. n ETA = 15, n ETA matched plasma = 14, n Sputum = 20, n Sputum matched plasma = 19, n BAL = 6, n BAL matched plasma = 3. d Normalized levels of cytokines/chemokines for COVID-19 respiratory and plasma samples separately. Red color indicates higher cytokine/chemokine levels. Source data are provided as a Source Data file. with FcγR binding activities were strongly featured in COVID-19 plasma compared to non-COVID-19 plasma (Fig. 4c).
Increasing cellular infiltrates in respiratory specimens during disease progression. To determine cellular immunity in respiratory specimens of COVID-19 patients, samples underwent multi-parameter flow cytometry and analysis using the Spectre R package 26 . Cells were clustered using Flow Self-Organizing Map (FlowSOM) 27 and plotted using Fast Interpolation-based t-distributed Stochastic Neighbour Embedding (FIt-SNE) 28 . Two flow cytometry panels were used to ensure accurate profiling of myeloid and lymphoid cell populations ( Supplementary Fig. 3a, 4a, b, Supplementary Table 7).
Clustering of respiratory samples in the myeloid panel revealed that CD66b + neutrophils dominated, with varying levels of CD16 expression (Fig. 5a). CD14 + macrophages and CD4 + and CD8 + T cells were also detected, but at lower frequencies. While the cellular component was variable across samples, CD16 hi and CD16 lo neutrophils were present in all COVID-19 patients apart from patient #043 (BMT recipient; Fig. 5a, Supplementary  Fig. 4c). Although there were only two COVID-19 patients (patients #026, #049) with multiple ETA samples, we still observed an increase in cellular infiltrates over time, including CD16 lo neutrophils (Fig. 5b). In the respiratory specimens of 6 non-COVID-19 patients, lower levels of neutrophils and macrophages were detected ( Supplementary Fig. 4c). Non-COVID-19 patient #059 had a large population of CD16neutrophils, while a high frequency of CD16 lo neutrophils was detected in blood, indicating a dominant immature neutrophil population in this patient ( Supplementary Fig. 4c, e).
After excluding neutrophils and monocytes/macrophages in respiratory samples, CD8 + T cells were the major population of lymphocytes, with varying levels of CD4 + T cells and natural killer (NK) cells (Fig. 5c, Supplementary Fig. 4d). Increasing infiltrates of T cells over time were found in patients #026 and #049 (Fig. 5d), similar to neutrophils. Interestingly, in patient #026, the lymphocyte population was dominated by NK cells early (d16 and d17) and T cells gradually infiltrated and dominated overtime (d23). Low lymphocyte levels were detected in fatal patient #021.
A volcano plot was generated to determine fold differences in immunological features between respiratory and blood samples (Fig. 5e, f). While cell numbers were higher in blood, higher frequencies of intermediate (CD14 + CD16 + ) monocytes/macrophages, activated (HLADR + CD38 + ) and EM-like (CD27 -CD45RA -) CD4 + and CD8 + T cells were found in COVID-19 respiratory specimens compared to blood. Respiratory specimens had a higher neutrophil to T cell ratio (Fig. 5f). Conversely, the ratio of CD4 + to CD8 + T cells was lower in respiratory samples (p = 0.0065), indicating a high prevalence of CD8 + T cells in respiratory specimens (Fig. 5f).
Overall, neutrophils (CD16 +/− ) dominated in the respiratory samples of COVID-19 patients, with varying levels of monocytes/ macrophages, T cells (CD4 + and CD8 + ), NK cells, and B cells. T cells in the respiratory samples exhibited an activated and EMlike phenotype compared to paired blood samples, with lower CD4 + to CD8 + T cell ratios.
IFN-α2 and IL-12p70 levels in ETA and RBD neutralizing activity in sputum negatively correlate with days of hospital stay. To understand associations between clinical features and serological responses in the respiratory specimens, correlations between clinical data (age, weight, height, BMI, days post disease onset, days of hospital stay) and serological features (cytokines and chemokines, sIL-6Rα, ADAMTS4, anti-RBD IgM, IgG, IgA and sVNT inhibition) were performed for ETA and sputum samples separately (Fig. 6a-d). IFN-α2 and IL-12p70 levels in ETA negatively correlated with days of hospital stay, albeit there were low levels of IL-12p70 in respiratory samples (Fig. 6a, b). In sputum, sVNT inhibition activity negatively correlated with days of hospital stay, but positively correlated with levels of MCP-1, IL-6, IL-8 and IL-10 ( Fig. 6c, d).
COVID-19 patients with higher NIH scores had more robust humoral immune responses in blood. While immune responses in blood samples between ward and ICU patients has been investigated in many studies, the classification of patients using NIH scores based on symptoms might correlate better with their immune responses. Unsurprisingly, more patients with higher NIH scores of 4-5 required ICU during hospitalization (Fig. 1d), while NIH scores of 2-3 were in the mild/moderate group. While all blood samples were collected during the acute phase of the infection, samples were grouped into hospital admission (V1) and hospital discharge (V7) for analyses. Although there were no differences in the overall cytokine/chemokine levels between the two NIH severity groups, IL-8 levels in the severe/critical group increased at V7 compared to V1 (p = 0.0004; Fig. 7a), indicating delayed or prolonged innate immune activation. Levels of sIL-6Rα were significantly higher in the severe/critical group than the mild/moderate group at both V1 (p = 0.027) and V7 (p = 0.0302), with the severe/critical group having higher sIL-6Rα and lower IL-6:sIL-6Rα ratios at V7 than V1 (Fig. 7a).
COVID-19 patients in the severe/critical group had comparable frequencies of immune cells, while they had lower T cell and eosinophil frequencies (p = 0.0011; p = 0.0473) than the mild/ moderate group at admission (V1; Fig. 7d). Interestingly, frequencies of mucosal-associated invariant T (MAIT) cells and ɣδ T cells negatively correlated with days in hospital (p = 0.0022 and p = 0.0024, respectively; Fig. 7d).
Overall, while cytokine levels were similar between the two severity groups, patients with more severe symptoms had more robust antibody responses towards the SARS-CoV-2.
Dexamethasone did not alter immune responses in COVID-19 patients in blood. Effects of dexamethasone, a corticosteroid anti-inflammatory drug, with/without remdesivir on immune responses in blood are unclear. We found very few differences in immune profiles between patients with/without dexamethasone. IL-8 and sIL-6Rα levels at discharge were significantly higher than at admission in the dexamethasone (with/without remdesivir) group, but similar levels were observed without treatment (Fig. 8a). Patients on treatment had lower anti-inflammatory IL-10 levels at discharge (p = 0.0281; Fig. 8a). Conversely, the humoral responses of patients receiving drugs were not compromised. Patients receiving dexamethasone (with/without remdesivir) generated robust SARS-CoV-2-specific antibody responses (Fig. 8b). Given that 29/33 severe/critical patients were on treatment, compared to 7/27 mild/moderate patients, high antibody levels in the drug group were likely due to disease severity rather than drug treatment. PLSDA revealed that patients prior to drug therapy had higher antibodies against the NP of human coronavirus OC43 rather than SARS-CoV-2, providing insights into potential drug treatment based on patient antibody responses at hospital admission (Fig. 8c). No significant differences were found in cellular responses, apart from lower T cell frequency in the drug group (Fig. 8d).

Discussion
Immunity to SARS-CoV-2 in the respiratory tract, the primary site of infection, is incompletely understood. We found differential inflammatory status in the respiratory tract and blood of COVID-19 patients, with high magnitude of MCP-1, IL-6, and IL-8 in respiratory specimens. While high SARS-CoV-2-specific IgG and IgM were detected in COVID-19 respiratory samples, IgG with FcγR-binding profiles were more prominent in blood. We found higher frequencies of neutrophils, intermediate CD14 + CD16 + monocytes, activated HLA-DR + CD38 + CD4 + T cells, EM-like CD4 + and CD8 + T cells in COVID-19 respiratory compared to blood samples. In blood, similar humoral immune responses were observed in patients with/without dexamethasone treatment.
High levels of cytokines are commonly found in the blood of COVID-19 patients [29][30][31][32] . In respiratory samples, variable cytokine levels (IL-10, IL-17A, IL-18) were detected, while monocyte chemoattractants (MCP-1, MIP-1α, MIP-1β) and innate cytokines (IL-6, IL-10) were at high levels 18,19 . We found hypercytokinemia in respiratory samples compared to blood, especially IL-6, IL-8, and MCP-1, indicating an inflammatory environment that attracts leukocytes, including neutrophils and monocytes 33,34 . Since most patients did not have similar cytokine profiles in blood and respiratory samples after normalizing the cytokine levels within each sample type, measuring both blood and respiratory inflammation, especially in sputum, might be needed to accurately determine the inflammatory status of the patients. Interestingly, IFN-α2 level in ETA samples negatively correlated with days of hospital stay in our cohort. Similarly, IFN levels in nasopharyngeal samples of COVID-19 patients negatively correlated with viral load 35 , indicating that viruses might be better controlled in the respiratory tract with higher IFN levels. Conversely, while higher plasma IL-12 levels were associated with more severe disease in COVID-19 patients 36 , ETA IL-12 levels also negatively correlated with days in the hospital, potentially due to its enhancement of CD8 + T cell activation 37 .
As an anti-inflammatory drug, dexamethasone can reduce proinflammatory cytokines including IL-6 and IL-8 44,45 . In blood, we found no significant difference in cytokine/chemokine levels between COVID-19 patients receiving dexamethasone and untreated patients. Although it has been speculated that dexamethasone can reduce the ability of B cells to produce antibodies 46 , we showed similar antibody levels in patients with and without dexamethasone. Therefore, severely-ill COVID-19 patients might benefit from dexamethasone treatment as reported [5][6][7] , and such treatment does not dampen humoral immunity.
There are limitations to the current study. Firstly, ETA samples were only collected from patients with severe disease requiring ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30088-y invasive oxygen support, therefore, it is unclear whether COVID-19 patients with milder symptoms had less robust immune responses in the respiratory tract. Additionally, most patients in the severe/critical group received dexamethasone, which could be an intercorrelating factor for the differences observed between severity groups. Moreover, while the non-COVID-19 controls provided insights onto the immune status in hospitalized individuals, the comparisons would benefit from larger numbers of non-COVID patients with more homogenous diseases.
Overall, innate and adaptive immune responses are generated in respiratory and blood samples of COVID-19 patients. While immunological features detected in the peripheral blood might predict clinical outcomes, monitoring immune responses in the respiratory samples can be of benefit prior to initiation of therapeutic interventions for COVID-19 patients given the disparity observed between respiratory and blood specimens.
Methods COVID-19 study participants and specimens. We enrolled 60 SARS-CoV-2 PCR-positive patients admitted to Austin Health (Victoria, Australia) and 6 PCRnegative patients as negative controls for SARS-CoV-2 antibodies. We enrolled an additional 24 COVID-19 patients through Austin Hospital, Royal Melbourne Hospital, Alfred Hospital and Westmead Hospital. Two COVID-19 patients and three SARS-CoV-2 PCR-negative patients died during the study. Participants of the current study were not compensated. Peripheral blood was collected in heparinized, ethylenediaminetetraacetic acid (EDTA) or serum tubes during hospitalization and centrifuged to collect plasma or serum. Peripheral blood monocular cells (PBMCs) were isolated via Ficoll-Paque separation. Single-cell suspensions were isolated from tissues as previously described 47,48 . ETA samples were obtained as part of routine suctioning of the endotracheal tube airway and involved the passage of a catheter for suctioning into a sterile respiratory sample trap. Sputum samples were spontaneously collected into a sterile container. Respiratory samples from the Alfred Hospital were residual samples taken as part of routine care. Pleural fluid was collected by thoracentesis as part of a routine procedure. The thoracentesis involved the percutaneous insertion of a catheter into the pleural space and collection of pleural fluid into a sterile container. Demographic, clinical and sampling information for COVID-19 patients are described in Supplementary  Table 1. For the 6 respiratory samples with undetectable cytokine/chemokine levels, we were still able to detect antibody levels, reflecting the high quality and integrity of the samples (Supplementary Table 6). Genomic sequencing and bioinformatic analysis. Extracted RNA from RT-PCR positive samples underwent tiled amplicon PCR and Illumina short-read sequencing, quality control, consensus sequence generation and alignment as previously described 49 . A single sequence per patient was used for phylogenetic analysis 22 , with a maximum-likelihood phylogenetic tree generated using IQ-Tree (v2.1.-, options "-mset GTR + G4 -bb 1000") 50 and visualized using the ggtree package (v.1.14.6) in R (v3.5.3) 51 . Genomic clusters were defined using a hierarchical clustering algorithm; genomic transmission networks grouped multiple clusters supported by epidemiological and genomic data.
Phenotypic whole blood immune analyses. Fresh whole blood (200 μl per stain) was used to measure CD4 + CXCR5 + ICOS + PD1 + follicular T cells (Tfh) and CD3 − CD19 + CD27 hi CD38 hi antibody-secreting B cell (ASC; plasmablast) populations as described 15,52 as well as activated HLA-DR + CD38 + CD8 + and HLA-DR + CD38 + CD4 + T cells, intermediate CD14 + CD16 + and classical CD14 + monocytes, activated CD3 − CD56 + NK cells, MAIT cells, ɣδ-T cells, as per the specific antibody panels (Supplementary Table 7; gating strategy is presented in Supplementary Fig. 3b, c). After whole blood was stained for 20 min at room temperature in the dark, samples were lysed with BD FACS Lysing solution (BD Biosciences, San Jose, California, USA), washed and fixed with 1% PFA. AccuCheck Counting Beads (Thermo Fisher Scientific, Carlsbad, CA, USA) were added for calculating absolute numbers just prior to acquisition. All samples were acquired on a LSRII Fortessa (BD) using the software BD FACS DIVA v8.0.1. Flow cytometry data were analyzed using FlowJo v10 software.  Table 7). After fixing with 1% PFA, the samples were acquired on a LSRII Fortessa (BD) using the software BD FACS DIVA v8.0.1. AccuCheck Counting Beads (Thermo Fisher Scientific) were added for calculating absolute numbers just prior to acquisition. Flow cytometry data were analyzed using FlowJo v10 software.
SARS-CoV-2 RBD ELISA. RBD-specific ELISA for detection of IgM, IgG and IgA antibodies was performed as previously described 31,53,54 , using flat bottom Nunc MaxiSorp 96-well plates (Thermo Fisher Scientific) for antigen coating (2 µg/ml), blocking with PBS with w/v 1% BSA and serial dilutions in PBS with v/v 0.05% Tween and w/v 0.5% BSA. Plates were read on a Multiskan plate reader (Labsystems, Vantaa, Finland) using the Thermo Ascent Software for Multiskan v2.4. Inter-and intra-experimental measurements were normalised using positive control plasma from a COVID-19 patient run on each plate. Endpoint titres were determined by interpolation from a sigmodial curve fit (all R-squared values >0.95; GraphPad Prism 9) as the reciprocal dilution of plasma that produced >15% (for IgA and IgG) or >30% (for IgM) absorbance of the positive control at a 1:31.6 (IgG and IgM) or 1:10 dilution (IgA). Seroconversion was defined when titres were above the mean titre (plus 2 standard deviations) of non-COVID-19 control respiratory or plasma samples.
Microneutralization assay. Microneutralization activity of serum samples was assessed as previously described 55 . SARS-CoV-2 isolate CoV/Australia/VIC01/ 2020 56 was propagated in Vero cells (ATCC #CCL-81) and stored at −80°C. Sera were heat-inactivated at 56°C for 30 min and serially diluted. Residual virus infectivity in the serum/virus mixtures was assessed in quadruplicate wells of Vero cells incubated in serum-free media containing 1 μg/ml of TPCK trypsin at 37°C and 5% CO 2 . The viral cytopathic effect was read on day 5. The neutralizing antibody titer was calculated using the Reed-Muench method 55 .

SARS-CoV-2 surrogate virus neutralisation test (sVNT)
. The plasma samples were tested in neat, and the respiratory samples were tested at 1:9 dilution or at their original dilutions for more diluted samples. The sVNT blocking ELISA assay (manufactured by GenScript, NJ, USA) was carried out essentially as described 54 , which detects circulating neutralizing SARS-CoV-2 antibodies that block the interaction between RBD and ACE2 on the cell surface receptor of the host. A HRP-conjugated recombinant SARS-CoV-2 RBD fragment bound to any circulating neutralizing anti-RBD antibodies preventing capture by the human ACE2 protein in the well, which was subsequently removed in the following wash step. Substrate reaction incubation time was 20 min at room temperature and results were read spectrophotometrically. Colour intensity was inversely dependent on the titre of anti-SARS-CoV-2 neutralizing antibodies.
Coupling of carboxylated beads. As previously described 25 , a custom multiplex bead array was designed and coupled with SARS-CoV-2 spike 1 (Sino Biological), spike 2 (ACRO Biosystems), RBD (BEI Resources) and nucleoprotein (ACRO Biosystems), as well as SARS and hCoV (229E, NL63, HKU1, OC43) spikes and Fig. 5 Higher frequencies of activated immune cells and EM-like CD4 + and CD8 + T cells in COVID-19 respiratory compared to paired blood samples. Flow Self-Organizing Map (FlowSOM) analyses of cellular content in the respiratory tract. a-d Metacluster of cells and expression level of markers in the a myeloid antibody panel and c lymphocyte antibody panel. Multiple endotracheal tube aspirate (ETA) samples from intensive care unit (ICU) patients #026 and #049 as well as one sputum sample from a fatal patient #021 were shown as example (b, d). e Volcano plot showing fold difference of 62 immunological features between paired respiratory and blood samples. f Comparisons of cellular immune features between respiratory and paired blood samples. n Respiratory = 14, n Respiratory matched blood = 13. Statistical significance was determined with a two-sided Mann-Whitney test. Dotted lines connect the most closely matched blood and respiratory samples from each patient. Colours indicate each patient. Source data are provided as a Source Data file.  Table 5). In addition, SARS-CoV-2 and HKU-1 spike trimers (kind gifts from Adam K. Wheatley), as well as SARS-CoV and NL63 spike trimers (BPS Bioscience) were also coupled. Tetanus toxoid (Sigma-Aldrich), influenza hemagglutinin (H1Cal2009; Sino Biological) and SIV gp120 (Sino Biological) were included in the assay as positive and negative controls respectively. Antigens were covalently coupled to magnetic carboxylated beads (Bio Rad) using a two-step carbodiimide reaction and blocked with 0.1% BSA, before being resuspended and stored in PBS 0.05% sodium azide.
Data normalization. For all multivariate analysis, Tetanus, H1Cal2009, and BSA antigens (positive controls) were removed, as well as SIV (negative control). Low signal features were removed when the 75 th percentile response for the feature was lower than the 75 th percentile of the BSA positive control. Right shifting was performed on each feature (detector-antigen pair) individually if it contained any negative values, by adding the minimum value for that feature back to all samples within that feature. Following this, all data were log-transformed using the following equation, where x is the right-shifted data and y is the right-shifted logtransformed data: y = log10(x + 1). This process transformed the majority of the features to having a normal distribution. In all the subsequent multivariate analyses, the data were further normalized by mean centering and variance scaling each feature using the z-score function in Matlab 2017b (MathWorks, Natick, MA). Plasma and respiratory samples were analysed separately. When analysing samples at time of hospital discharge, to adjust for the confounder of time from symptom onset, each of the features were iteratively regressed with ordinary least squares regression, using the residuals as input for the analysis 57 .
Feature selection using elastic Net/PLSDA. To determine the minimal set of features (signatures) needed to predict categorical outcomes (COVID-19 diagnosis, NIH scores, drug therapies), a three-step process was developed 58 . First, the data were randomly sampled without replacement to generate 2000 subsets. The resampled subsets spanned 80% of the original sample size, or sampled all classes at the size of the smallest class for categorical outcomes, which corrected for any potential effects of class size imbalances during regularization. Elastic-Net regularization was then applied to each of the 2000 resampled subsets to reduce and select features most associated with the outcome variables. The Elastic-Net hyperparameter, α, was set to have equal weights between the L1 norm and L2 norm associated with the penalty function for least absolute shrinkage and selection (LASSO) and ridge regression, respectively 59 . By using both penalties, Elastic-Net provides sparsity and promotes group selection. The frequency at which each feature was selected across the 2000 iterations was used to determine the signatures by using a sequential step-forward algorithm that iteratively added a single feature into the PLSDA model starting with the feature that had the highest frequency of selection, to the lowest frequency of selection. Model prediction performance was assessed at each step and evaluated by 10-fold cross-validation classification error for categorical outcomes. The model with the lowest classification error within a 0.01 difference between the minimum classification error was selected as the minimum signature. If multiple models fell within this range, the one with the least number of features was selected and if there was a large disparity between calibration and cross-validation error (over-fitting), the model with the least disparity and best performance was selected.
PLSDA. Partial least squares discriminant analysis (PLSDA), performed in Eigenvectors PLS toolbox 8.2 in Matlab 2017b, was used in conjunction with Elastic-Net, described above, to identify and visualize signatures that distinguish categorical outcomes (COVID-19 diagnosis, NIH scores, drug therapies). This supervised method assigns a loading to each feature within a given signature and identifies the linear combination of loadings (a latent variable, LV) that best separates the categorical groups. A feature with a high loading magnitude indicates greater importance for separating the groups from one another. Each sample is then scored and plotted using their individual response measurements expressed through the LVs. The scores and loadings can then be cross-referenced to determine which features are loaded in association with which categorical groups (positively loaded features are higher in positively scoring groups, etc.). All models go through 10-fold cross-validation, where iteratively 10% of the data is left out as the test set, and the rest is used to train the model. Model performance is measured through calibration error (average error in the training set) as well as crossvalidation error (average error in the test set), with values near 0 being best. All models were orthogonalized to enable clear visualization of results. PLSDA scores and loadings plots were plotted in Prism v8.
Hierarchical clustering. We visualized the clustering of DRASTIC respiratory and blood samples based on only SARS-CoV-2 antigens or all features using unsupervised average linkage hierarchical clustering of normalized data using MATLAB 2017b. Euclidean distance was used as the distance metric. sIL-6Rα and ADAMTS4 ELISAs. Soluble protein levels were all measured using DuoSet ELISA kits for each protein (R&D Systems, Minneapolis, MN, USA) according to the manufacturer's instructions. DuoSet ELISA ancillary reagent kit (R&D Systems) was used for respiratory fluids and in-house reagents with the same composition were used for plasma samples. In brief, 96-well R&D ELISA microplates (respiratory fluids) or 96-well Nunc Maxisorp ELISA plates (ThermoFisher, plasma) were coated with capture antibody overnight, followed by blocking with 1% w/v BSA for a minimum of 1 h. Samples and standard proteins were added and incubated for 2 h at room temperature, followed by detection antibody for a further 2 h. Lastly, streptavidin-HRP, substrate solution and stop solution (2 N H 2 SO 4 ) were added subsequently for 20 min each. Plates were read on a Multiskan plate reader (Labsystems) using the Thermo Ascent Software for Multiskan v2.4. Plasma samples were diluted in 1:300 for sIL-6Rα ELISAs. Respiratory fluids were diluted in 1:50/1:150 for sIL-6Rα ELISA accounting in the original dilution factors and tested without further dilution for ADAMTS4 ELISA.
Computational flow cytometry analysis. Computational analysis of data was performed using the Spectre R package (v0.4.1) 26 (https://github.com/ImmuneDynamics/ Spectre). Samples were initially prepared in FlowJo, and populations of interest were exported as CSV files containing raw (scale value) data. In R (v4.0.2), data were subject to arcsinh transformation and clustering using FlowSOM (v1.20.0) 27 . For visualisation, cells in the myeloid panel were subjected to sample-weighted downsampling based on absolute cells/uL counts in the blood, whereas cells in the lymphoid panel were unchanged to preserve samples with low cell numbers. Cells from the lymphoid panel, and downsampled cells from the myeloid panel, were then distributed in 2D via dimensionality reduction (DR) using Fast Interpolation-based t-distributed Stochastic Neighbour Embedding (FIt-SNE, v1.2.1) 28 . Data from both the lymphoid and myeloid panel were subject to two rounds of clustering and DR. The initial round of clustering and DR was used to filter out cellular debris and non-immune cells exhibiting high autofluorescence, using arcsinh transformed expression of CD45, CD3, CD4, CD8, CD14, CD16, CD19, CD64, CD66b, CD38, HLA-DR, and Live-Dead for the myeloid Fig. 7 COVID-19 patients with higher NIH scores display more robust humoral immune responses towards SARS-CoV-2 in blood. a Levels of cytokines, soluble IL-6 receptor α (sIL-6Rα), and IL-6:sIL-6Rα ratio, b anti-RBD IgM, IgG, and IgA titres, microneutralization titres and c Partial Least-Squares Discriminant Analysis (PLSDA) scores and loadings plot of antibodies against human coronavirus between mild/moderate and severe/critical COVID-19 patients. d Volcano plot showing fold difference of 83 immunological features in blood samples between mild/moderate and severe/critical COVID-19 patients, and comparisons of cellular subset frequencies and correlation with days stayed in hospital. n Mild-Moderate V1 = 25, n Mild-Moderate V7 = 14, n Severe-Critical V1 = 31, n Severe-Critical V7 = 16. The bounds of the box plot indicate the 25th and 75th percentiles, the bar indicates medians, and the whiskers indicate minima and maxima. Statistical significance was determined with a two-sided Kruskal-Wallis test followed by Dunn's multiple comparisons test. Partial Least-Squares Discriminant Analysis was performed for antibodies measured with multiplex bead array assay. Volcano plots were created using a two-sided Wilcoxon rank-sum test and statistics were corrected with FDR adjustment. Correlation was determined with Spearman's correlation. V1, hospital admission; V7, hospital discharge. Source data are provided as a Source Data file.
Volcano plots and heatmaps were created using the Spectre R package 26 , where comparisons were performed using a Wilcoxon rank-sum test (equivalent to the Mann-Whitney test) with the wilcox.test function in R. Statistics displayed in volcano plots were corrected with a False Discovery Rate (FDR) adjustment.
Statistical analyses. Statistical significance was assessed using two-sided Mann-Whitney, Wilcoxon signed-rank test or Kruskal-Wallis test with Dunn's correction for multiple comparisons in Prism 9 (GraphPad) unless stated otherwise. Correlations were assessed using two-tailed Spearman's correlation coefficient (r s ) and visualized in R v3.6.2 as heatmaps using the corrplot package or using the online Morpheus heatmap software (https://software.broadinstitute.org/morpheus; the Broad Institute, MA, USA) and p-values of correlations were corrected for multiple comparisons by FDR in R v3.6.2. P-values lower than 0.05 were considered statistically significant.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All data generated or analysed during this study are included in this published article (and its supplementary information files). A Source Data file is provided with this paper. All relevant data are also available from the authors. The viral sequences isolated from nasal swabs that support the findings of this study are available on the GISAID database with ID numbers provided in the Source data file. Registration to access the database is free and open to anyone using the link: https://www.gisaid.org/registration/register/.